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We discuss the possibility of sampling exponential moments 
of the canonical phase from the s-parametrized phase space 
functions. We show that the sampling kernels exist and are 
well-behaved for any s > — 1, whereas for s = — 1 the kernels 
diverge in the origin. In spite of that we show that the phase 
space moments can be sampled with any predefined accu- 
racy from the Q-function measured in the double-homodyne 
scheme with perfect detectors. We discuss the effect of imper- 
fect detection and address sampling schemes using other mea- 
surable phase-space functions. Finally, we discuss the prob- 
lem of sampling the canonical phase distribution itself. 

PACS number(s): 42.50.Dv 



I. INTRODUCTION 

Studying the role of phase in quantum mechanics has a 
long history (for a review on the phase concepts, see, e.g., 
jij). Its importance in today's problems is also appar- 
ent. For example, phase in atomic systems has recently 
been used for storing quantum information , and phase 
and photon number measurements have been considered 
as a basis in some quantum teleportation schemes ||. 
Notwithstanding the various phase-dependent effects in 
quantum physics, phase itself has not been uniquely mea- 
sured and its very definition as physical quantity has been 
subject to many disputes. Whereas for highly excited 
(quasi-classical) states different approaches give similar 
results, the various concepts differ in the phase properties 
of quantum states close to vacuum. Therefore the ques- 
tion has been arisen of what are the differences between 
these approaches and how relevant are they experimen- 
tally. In this paper we concentrate on the canonical phase 
and its relation to s-parametrized phase-space functions, 
with special emphasis on the measurability of its expo- 
nential moments by "weighted" averaging of measured 
phase space functions. 

The canonical phase distribution P(ip) of a radiation 
field mode (harmonic oscillator) prepared in a quantum 
state described by a density operator g is defined by 



(2) 
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P{<p) = (2ir)- 1 (ip\e\^), 
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Even though there has been no known experimental 
scheme that is directly governed by P(tp), this distribu- 
tion has very nice properties: it is non-negative, conju- 
gated to the photon-number distribution (in the sense 
that a phase shifter shifts a phase distribution while 
a number shifter does not change it there exist 

number-phase uncertainty relations |^|, and in compari- 
son to other phase distributions, P(<p) is the most sharp 
one. 

The lack of direct experimental availability of P(<p) has 
led us to the search of schemes for sampling the canoni- 
cal phase statistics from quantities that can be measured 
directly |Hq]. In balanced homodyne detection (for a 
review on quantum state measurement using homodyn- 
ing, see, e.g., ||), the exponential moments ^ of the 
canonical phase, 
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where the Fock state expansion of the (unnormalizable) 
phase states \ip) reads 



can be sampled by integrating the measured quadrature- 
component statistics multiplied by well-behaved kernel 
functions |6]-[|]. An advantage of the method is that it 
applies to both the quantum regime and the classical 
regime in a unified way. Of course, the question has been 
as of whether or not it is possible to find other (and pos- 
sibly better) measurement schemes suitable for sampling 
the exponential moments of the canonical phase. 

It is well known that balanced double-homodyne de- 
tection (eight-port homodyning, p0|) provides us with a 
two dimensional set of data whose statistics correspond 
to a s-parametrized phase-space function W s (q,p) with 
s < — 1 |n|. In this scheme, the limiting case of s = — I , 
which corresponds to the Husimi Q-function Q(q,p) = 
W-i(q,p), requires perfect detection, i.e., 100% detec- 
tion efficiency. Having a sampling scheme leading from 
a measured s-parametrized phase-space function to the 
exponential canonical-phase moments would be the most 
direct method of measuring the exponential moments of 
the canonical phase. Since each measurement event (q,p) 
already yields a phase value Arg(g + ip), the measured 
values only need to be "weighted" by the kernel func- 
tions in the averaging procedure yielding the exponential 
moments ^k- 

There are also measuring schemes, e.g., unbalanced 
homodyning, suitable for determining s-parametrized 
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phase-space functions W s (q,p) with larger values of s 
fjjfl . However, in these schemes the functions W s (q,p) 
are not obtained in terms of the statistics of measure- 
ment events (q,p), but they are obtained pointwise for 
each phase-space point (q,p) set up in the experiment. 
Moreover, they are typically reconstructed from the mea- 
sured data rather than measured directly. Nevertheless, 
it is interesting to ask the question of the prospects of 
phase measurement in schemes of that type. 

In this paper, we try to answer the questions raised 
above, focusing our attention to the problem of using bal- 
anced double-homodyne detection for sampling the ex- 
ponential moments of the canonical phase. In Sec. || we 
present the kernels that relate the s-parametrized phase- 
space functions to the exponential phase moments, and in 
Sec. Ill we apply the results to direct sampling of the ex- 



ponential phase moments in balanced double-homodyne 
detection. Other measurement schemes are discussed in 
Sec. IV. Section [v| addresses the problem of determining 
the phase distribution itself, and a conclusion is given in 



Sec. VI 



II. THE KERNEL FUNCTION 

Our task is to find the kernel function Kk(r;s) such 
that the exponential moments of the canonical phase can 
be given by (k > 0) 



* fe = (E k ) 
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dpe lkip / rdrW(r,p;s)K k (r;s), (4) 
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and = \&£, where 



E=J2\n){n+l\ 



(5) 



n=0 



In Eq. (0), the phase-space function W(r, p; s) is written 
in polar coordinates, i.e., W(r, ip; s) = W s (r cos ip, r simp). 
Note that e lklp K k (r; s) is the (— s)-parametrized phase- 
space function of the operator E k . We now take advan- 
tage of the expression jjj 
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where the expectation value of the normally ordered cor- 
relations of the photon creation and destruction opera- 
tors can be calculated by means of W(r, ip; s) as [JJ 
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(L k , Laguerre polynomial). Combining Eqs. (^), and 
(0) , we derive (Appendix 0) 
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Here, the function fl^ k \p 2 ) is given by 

rt k Xp 2 ) = e-P 2 Y. t T^p^ 1 (9) 
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where 



4 fe ) = / dip, sin fe - 2 p x 



...J dp>i sin fc 1 1 p l ...J dipk-i | [sin 2 ip x 

x (l + sin 2 ^ 2 (l + ...(l + sin 2 ^-i)))] n }. (10) 

It is worth noting that Kk(r;s) is unique, which follows 
from the fact that Kk(r;s) is the phase-space function 
of E k and from the uniqueness of phase-space represen- 
tations. This is in contrast to the kernel functions that 
relate quantities to the quadrature-component statistics 
measured in balanced homodyne scheme, where certain 
functions can be added to the kernels without changing 
the result 010. 

The integral in Eq. (@) converges for s > — 1 because 



|ft (fe V)| <e-"X 



(11) 



14 being some constant. Plots of the kernel function for 
different values of s and k are shown in Fig. |l|. We can see 
that Kk(r; s) monotonically increases with r from zero to 
one for s > 0. If s < 0, then Kk(r; s) attains the maximum 
at a finite value of r. The position of the maximum shifts 
towards the origin and the value of the maximum tends 
to infinity as s— > —1. Hence the kernels that relate the 
exponential phase moments to the Q-function diverge at 
r = 0. To be more specific, it can be shown (Appendix 
B) that 



K k (r; -1) cx r 



-k 
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near the origin. 

Though the function Kk{r) = K).(r;—1) diverges, it 
can be used to obtain the exponential phase moments 
from the Q-function Q(r, p) = W(r, tp, —1). It is not 
difficult to prove that Eq. (0) can be rewritten as 



rdrQ k (r)K k {r) 



(13) 
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FIG. 1. The kernel function K k (r;s) for k = 1,2,3 and 
s = 0.75 (full line), s = (broken line), s = -0.75 (dash-dotted 
line), s = — 1 (dotted line). 

where 
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(Pn+fe,n = (n+k\g\n)). It follows from Eq. (|l4|) that Qk(r) 
(xr k for small r, and thus Q&(r) exactly compensates for 
the divergence of Kk(r), Eq. (|l2|). In other words, if 
the Q-function of the state is known exactly, then the 
integration in ( |l3| ) can be performed straightforwardly, 
thus yielding the sought However, measurement of 
the Q-function is always associated with some error, so 
that the region close to the origin of the phase space 
needs careful consideration in praxis. 
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FIG. 2. Double-homodyne scheme [^). The signal beam 
is split on a beam splitter BSi and the resulting beams are 
mixed with strong coherent local oscillator (LO) on the beam 
splitters BS2 and BS4. The LO beams at BS2 and BS4 stem 
from a common source, split at BS3, and their phases differ 
by 7r/2, determined by the A/4 phase shifter. The difference 
of photocurrents measured at the detectors Di and D2 is pro- 
portional to q and the photocurrent difference at D3 and D4 
is proportional to p. 



III. CANONICAL PHASE FROM DOUBLE 
HOMODYNING 

A. Statistical error 

Let us consider balanced double-homodyne detection 
(Fig. ||) and first assume perfect detection. Each exper- 
imental event then gives a pair of real numbers which, 
after rescaling, define a point in the phase space of the 
signal, and the probability density of detecting the space 
points is equal to the Q-function of the signal state [jll) . 
When the jth measurement yields the phase-space point 
with polar coordinates [vj,(pj) and altogether TV mea- 
surements are performed, then the exponential phase mo- 
ments can be estimated to be 
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— \ ]exp(k<pj)K k (rj 
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In order to answer the question of how close is to 
the actual moment ^Pfc, we calculate the mean value and 
dispersion of the estimate ( [j~5| ) over all possible measure- 
ment results. Since individual measurement outcomes 
are independent of each other, we can take advantage of 
the summation rule for mean values and dispersions of 
independent quantities. Thus, for the real part of x I , ^. est ^ 
we get the mean value 

1 N 

3=1 
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iVRe* fe -Re* fe , (16) 
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N 

as it should be, and the dispersion 

, N 

D(Re*i ost) ) = — ^[cos^)^- 
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= -^<j J dip rdr cos 2 (kip)Kl(r)Q(r,ip) 



dip / rdr cos(kip)K k (r)Q(r,ip) 
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Similar expressions hold for the imaginary part of \I>jj. 
Since cos 2 (fc<p) = 1/2 + cos(2kip) /2, after performing the 
angular integration in the first term on the right-hand 
side of Eq. (|l7j), the radial part contains the product 
Qo{r)K k (r) so that this integral over the divergent ker- 
nel can become infinite. Let no be the number of photons 
at which the Fock expansion of the state starts. Taking 
into account Eq. (|l4|), we see that the integrand behaves 
as oc r 2 ("o-'0+i f or small r. Thus, the exponential phase 
moments ^ k can be directly sampled from the double- 
homodyne data, provided that k < no + 1, because in this 
case the dispersion of the estimation is bounded. In the 
opposite case of k > no + 1, the statistical fluctuation di- 
verges so that the exponential phase moments cannot be 
sampled without a proper regularization of the kernels. 
Note that for states that contain the vacuum, regulariza- 
tion of the kernels is necessary for all exponential phase 
moments. 



B. Kernel regularization and sampling algorithm 

Since the main part of the statistical error arises from 
data close to the origin, it is natural to modify the pro- 
cedure by omitting the data falling inside a small cir- 
cle r<ro (see Fig. |^). Of course, such a deliberate 
data filtering introduces into the measurement a state- 
dependent systematic error. Nevertheless, the statistical 
error is reduced and the total error may be acceptable. 
Replacing Kk(r) by the regularized function K' k (r) ac- 
cording to 



K' k (r) = 6(r-r )K k (r) 



(18) 



[6{x), Heaviside step function], the systematic error of 
the fcth moment can be given by 



dipe ikv / rdrQ(r,p)K k (r) 
ii Jo 

ro 

rdrQ k (r)K k {r). (19) 
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FIG. 3. Output of simulated double-homodyne detection 
of a coherent state \a), a — 1 (a), and enlarged detail of the 
output around the origin (b). 

A measure of the total error is then the sum of the sta- 
tistical and systematic errors, 



Re a 



(tot) 



Re4 sys) 



D(Re *i ost) ) 



1/2 



(20) 



and Im a k tot ^ accordingly. 

From the example in Fig. [| it is seen that the statisti- 
cal error decreases with increasing radius ro [Fig. |](a)], 
whereas the systematic error increases with the radius 
[Fig. [|(b)]. The total error has thus a minimum at a 
certain radius [Fig. 0(c)], which can be regarded as the 
optimal radius for regularization. Unfortunately, the de- 
termination of the systematic error requires knowledge of 
the state. Nevertheless, an upper bound of the system- 
atic error can be estimated, without any a priori knowl- 
edge of the measured state. Assuming r < 1, we may 
write 



\Qk(r)\ =2£ 
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FIG. 4. Statistical (a), systematic (b), and total (c) errors 
of the real part of the sampled exponential phase moment 
Re^2 of a coherent state \a), a = 1 for different numbers of 
recorded events N. 



where we have used the inequality 
implied by positive definiteness of p. Hence, an upper 
bound of the systematic error can be estimated. Using 
Eq. ©, we find that 



r ( s y s ) i 



< 



dr r 



fc+i 



K k (r). 
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A typical state for which | cr£ SyS ' > | is of the order of magni- 
tude of upper-bound value is \ipk) = (|0) + \ k))/y/2. For 



this state, Qfe(r) = r k exp(— r 2 )/\ 



which yields one 



half of the upper bound value. Taking into account that 
Kk{r) cx l/r fe for r -C 1, we find from the inequality ( p^ ) 
that the upper bound of |crj?, sys ''| increases quadratically 
with r . The dependence on r of the upper bound of 

the | cr£ fays ' > | is shown in Fig. || Notice that the systematic 
error is smaller for higher k. 

The state-independent upper bound of the systematic 
error and the estimated statistical error can now be used 
to determine the upper bound of the total error. Its 
minimum then determines an appropriate regularization 
radius ro- A possible algorithm for optimized data pro- 
cessing is the following one. In the zeroth step, sampling 
of the desired exponential phase moments from all N 
measurement events is performed. Since also data with 
very small r may contribute to the result, the statistical 
error can be very large. In contrast to standard sam- 
pling technique, where there is no need for data storage, 



0.5 



systematic error 




FIG. 5. Upper bound of jo-^. sys '| estimated from the in- 
equality ( p^ ) for k— 1 (solid line), k = 2 (dashed line), fc = 3 
(dot-dashed line), and fc = 4 (dotted line). 



here the data within a certain small circle are stored. 
The radius of the circle should be slightly larger than the 
expected regularization radius. The regularized kernel 
function (18) is now used, with r being increased step 
by step, so that in the nth step n events closest to the 
origin are covered by ro . In each step statistical and sys- 
tematic errors are estimated. The value of ro for which 
the total error is minimized is used for calculation of the 
final result. 

Let us mention that the detrimental effect of diver- 
gent kernels K k (r) at r = (in connection with nonzero 
Q-function) resembles the experiment in |1C|] , where the 
statistics of sine and cosine phases are obtained from low- 
efficiency double homodyning. In the experiment, data 
giving rise to divergences are disregarded, which is criti- 
cized in p5[ from the argument that the disregarded data 
represent an extra noise in the statistics. In our case, 
we disregard data leading to high statistical error and 
include the resulting systematic error into the sampling 
scheme. 



C. Total error and number of measurements 

Let us assume that a particular phase moment "Ffc is 
desired to be determined with a prescribed total preci- 
sion cr!? • What is the necessary number of measure- 
ment events N? If there were no need for regularization 
and the precision were limited only by (finite) statistical 

fluctuation, then iVcx (c[ t ° t ' > ) _2 . When the vacuum con- 
tributes to the state to be measured and a regularization 
radius ro is introduced, then the total error reads 

at t) =A 1 (-lnr ) 1 / 2 N- 1 / 2 + B 1 r 2 , (23) 
at° t] = A k rl k N-^ + B k r 2 , k > 2, (24) 

where A k and B k are constants. The optimal regulariza- 
tion radius rQ° pt ^ , which minimizes the total error ( [24] ) 
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FIG. 6. Sampled exponential phase moments of a coherent 
state \a), a = l; (a) real part of ^jf 8 *', (b) imaginary part of 
*^ cst) . The error bars indicate the estimated statistical error. 
In the computer simulation, N — 10 6 events are recorded and 
perfect detection is assumed. The dashed regions correspond 
to the phase moments of the radially integrated Q-function. 



depends on N as 

r<° pt) oc JV-VPd+k)] . 
From this expression and Eq. ( p4| ) we find that 

^(tot) K N -l/(l+k) 



N oc 



(tot; 



-(i+fe) 



(25) 



(26) 



(27) 



(fc > 2). The case k = 1 needs separate consideration, be- 
cause of the logarithm, which does not provide us with 
a simple analytical expression. Obviously, N increases 
faster than (<7 1 tot ' ) ) -2 with decreasing error. Thus, we can 
see that in the limit of small total error ordinary homo- 
dyning (which does not require regularization) is better 
suitable for sampling exponential phase moments than 
the double homodyning, because it requires a smaller 
amount of data to achieve the same precision. 
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Re *j; est) 


Im*[ eBt) 


' u 


1 


0.7732 


0.7790 ± 0.0006 


-0.0008 ± 0.0007 


0.007 


2 


0.4805 


0.483 ±0.003 


0.001 ± 0.003 


0.061 


3 


0.2559 


0.26 ±0.01 


0.00 ±0.01 


0.160 


4 


0.1209 


0.13 ±0.02 


0.01 ±0.02 


0.277 



TABLE I. Comparison of the sampled exponential phase 
moments shown in Fig. |^ with the exact ones. The displayed 
optimized regularization radii ro refer to Re values cor- 
responding to Im vffe are similar in magnitude. 



k 




Re*< cst) 


Im *jf Bt) 


1 


0.7732 


0.7736 ± 0.0004 


0.0002 ± 0.0006 


2 


0.4805 


0.4795 ± 0.0009 


-0.0003 ±0.001 


3 


0.2559 


0.2573 ±0.0017 


0.0005 ± 0.0017 


4 


0.1209 


0.1200 ±0.0021 


0.0002 ± 0.0021 



TABLE II. Comparison of the exponential phase moments 
of a coherent state \a), a = 1, sampled in homodyne detec- 
tion with the exact ones. In the computer simulation, the 
27r-phase interval of the quadrature components is divided 



into 120 equidistant values, and altogether N - 
recorded. 



10 events are 



difference between the concept of canonical phase and 
the phase concept based on the radially integrated Q- 
function. 

In order to compare double homodyning with ordinary 
homodyning, we have also simulated homodyne detection 
of the quadrature-component statistics for sampling the 
exponential phase moments of the same coherent state as 
in the simulated double- homodyne experiment, using the 
method in H ||. The results are presented in Table |J. 
Comparing Tables [j] and [fj|, we see that (f° r equal to- 
tal numbers of events) the error in ordinary homodyning 
is indeed smaller than in double homodyning. The dif- 
ference between the errors observed in the two schemes 
increases with increasing index k of the moment. Note 
that for k = 4 the error in the double-homodyne scheme 
is ten times larger than in the ordinary homodyne mea- 
surement. 



E. Imperfect photodetection 



D. Computer simulation 

To demonstrate the feasibility of the method, we have 
performed Monte Carlo simulations of double-homodyne 
detection of the Q-function for sampling the exponential 
phase moments of a coherent state. Results are shown 
in Fig. ^ and Table |. From Fig. ^| and Table | it is seen 
that the sampled exponential phase moments are in good 
agreement with the exact ones. Note the strong increase 
of the error with the index k of the moment (for a detailed 
error analysis, see Fig. Further, a comparison between 
the dashed and undashed bars in Fig. H clearly shows the 



In a real experiment, the overall detection efficiency r] 
would be always smaller than 100%, but it can be very 
high, e.g., 77 = 99%. Nonperfect detection introduces ad- 
ditional noise into the sampling scheme and gives rise to 
an additional systematic error, which cannot be dimin- 
ished by increasing the number of measurements. The ef- 
fect of nonperfect detection is that the exponential phase 
moments of a "smoothed" quantum state are sampled 
rather than those of the true one. Since the additional 
noise is Gaussian, the phase-space function that is ac- 
tually recorded is not the Q-function but the function 
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Wi-2))- 1 (q,p). The sampling of ^ k from this function 
by means of the kernel function K k (r; —1) is equivalent 
to sampling of ^> k from the Q-function by means of the 
kernel function K k [r; — 1 + 2(?7 _1 — 1)]. For a given quan- 
tum state, the systematic error can thus be given by 

/>OC 

A r) * fc = / rdrQ k (r)[K k (r;-l)-K k (r;-3 + 2r 1 - 1 )]. 
Jo 

(28) 

Its result is the underestimation of the magnitude of the 
moment. Since the difference of the kernel functions is 
essentially nonzero only around the origin r = 0, the sys- 
tematic error A v ^/ k will be highest for states close to vac- 
uum. After a proper regularization, one can use Eq. ( psj ) 
to get a reasonable estimation of the error by substituting 
the measured statistics for the unknown Q-function. 
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IV. OTHER MEASUREMENTS OF THE 
PHASE-SPACE FUNCTIONS 

The phase moments ^ k can be also obtained from qua- 
sidistributions reconstructed in unbalanced homodyning 
p2[ , or cavity measurements jl6|,|l7]]. However, these 
methods do not yield W s (q, p) as statistics of events (q, p) , 
but the functions W s (q,p) are determined pointwise. The 
restriction in practice to a selected finite number of points 
necessarily results in a systematic error, because the in- 
tegration over the phase space is replaced by summa- 
tion over a finite number of points selected by the exper- 
imentalist. Having determined W s (q,p), the exponen- 
tial phase moments ^ k can then be reconstructed from 
W s (q,p) on the basis of Eq. (Q). Since the kernel func- 
tion K k (r; s) is well behaved for s>— 1, no problems with 
divergences arise here. 

In unbalanced homodyning displaced Fock-state distri- 
butions p(n,a) are measured p"2| , |l8[ . It can be expected 
that the cumbersome way of reconstructing from 
p(n, a) via W s (q,p) may be avoided and the reconstruc- 
tion can be performed directly from the measured data. 
This could be done in a similar way as in the reconstruc- 
tion of the density matrix in the Fock basis p9| . A sim- 
ilar approach can be used for different physical systems: 
statistics of the displaced Fock states of vibrating trap 
ions has been obtained in state-reconstruction experi- 
ments pQ| ) and schemes based on displaced Fock statis- 
tics of the cavity fields have been suggested Jl^,|l^]. In 
particular, the scheme of |l6| directly yields the Wigner 
function, from which the exponential phase moments can 
be obtained in a very straightforward way. Even though 
many interesting problems are related to these schemes, 
we will not deal with them in this paper in any more 
detail. 



FIG. 7. Plot of the function F(r, </j; 0.5). 



V. PHASE DISTRIBUTION 

In order to answer the question of the possibility of 
direct sampling of the phase distribution itself, we have 
first to answer the question of the existence of kernels 
F(r, ip — ip; s) such that 

/>27T />00 

P((p)= dip rdrW(r 7 iP;s)F(r,ip-iP;s). (29) 
Jo Jo 

Obviously, F(r, ip — ip; s) is the (— s)-parametrized phase- 
space function of the phase state \ip) in Eq. (||). In pl| ] 
it is shown that this function can be given by 



F(r,ip;s) = J2B m Ar,s)e^ m ' n ^ , 



(30) 



where 



B mn (r, s) 




m+l 



XL 



(31) 



for m>n, and B nm = B mn . The series (|30j) only converges 
for s > (for the limiting case s = 0, see also p^|). 

Let us express F(r, tp; s) for s > in terms of K(r; s). 
From Eq. (H) it follows that 



P 



^ OO 

(?) = ^ £ *« 



-ikep 



(32) 



fe=-c 



Combining Eqs. (|32]) and (|4|) and recalling Eq. (|29|), we 
may write F(r, ip — tp;s) in the form 



1 

F{r, p; s) = — 



l + 2^2K k (r;s)cos(kip) 



k=l 



(33) 
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When r — > oo then K k (r;s) — > 1, and thus F(r, ( p, s) 
— > <5(<^). Regrouping the terms in Eqs. (Al) and ( |A2| ) 
(for s > 0) and using a summation formula for Laguerre 
polynomials p3], we can rewrite K k (r; s) as 



Kfc(r;s) = r* 



s+1 
(-1)" 



fc+i 



,-2r 2 /(l+s) 



^ V(n+l)...(n + fc) Vl + 



4r" 
1^ 



, (34) 



which is suitable for computing F(r,tp;s). An example 
is displayed in Fig. [?]. 

The fact that for a large class of states W(r, <p; s) does 
not exist as a regular function for s > limits the ap- 
plicability of Eq. (|29|). Nevertheless, there exist states 
for which W(r, ip; s) for s > is a regular function which 
can be sampled using unbalanced homodyning. However, 
for s > the statistical error of the sampled Wjr, ip; s) in- 
creases with r , so that application of Eq. ( p9|) requires 
special regularization. 



dyning where a compensation of imperfect detection is 
possible for efficiencies down to ry > 0.5. (vii) Thus, in re- 
ply to the question posed in the Introduction, it does not 
seem that phase-space measurements based on doublc- 
homodyning are closer to canonical-phase measurement 
than quadrature-component measurements based on or- 
dinary homodyning. (viii) In contrast to ordinary ho- 
modyning however, the sampling functions in double ho- 
modyning are uniquely defined. This follows from the 
fact that they are actually phase-space representations 
of quantum mechanical operators, (ix) The exponen- 
tial phase moments can also be inferred from the data 
recorded in other schemes such as unbalanced homodyn- 
ing, in which s-parametrized phase-space functions are 
reconstructed pointwise. (x) Kernels for sampling the 
distribution of the canonical phase exist as regular func- 
tions only for s > 0. Even though for some states the cor- 
responding phase space functions exist and can be mea- 
sured using unbalanced homodyning, the behavior of the 
statistical error would require a special regularization of 
the scheme to be applicable. 



VI. SUMMARY AND CONCLUSION 
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APPENDIX A: SAMPLING KERNELS AND THE 
FUNCTION SI 

Substituting Eq. (0) into Eq. (||) and comparing with 
Eq. (Eh, we can express K k (r; s) as 



1=0 

where the coefficients 
I 



1-sY „ {k)rk f 2r 2 



1- s 



-ii 



n=0 V( n + l )---( n + k ) K " 

can be rewritten as 



(Al) 



(A2) 



with 



(fe) 



1 



T k/2 



dt k e- kt *z k , (A3) 



z k = l-e-ti, p| = 



(A4) 



Substituting this expression into Eq. ( |A1[ ) and using the 
summation rule 
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1=0 



= [l-Zk- 



1 - s 

fc-i 



exp 



z fc (l-«)/2-l 



(A5) 



we arrive at 



K k (r;s) = — WT \ dhie^... dt k e 

2(l-e-pl)r 2 



l + .s+(l-.s)e" p '-- 



-fc-i 



exp 



l + s + (l-s)e~ p fc 



(A6) 



Note that the series in Eq. (A5) is only convergent for 
\z k (l - s)/2\ < 1. We have z k < 1, thus 



(1 — s)/2| < 1 => s>-l 



(A7) 



must hold so that the Q-function (,s = — 1) represents 
limiting case for sampling the phase moments. 

The multiple integration in Eq. ( |A6| ) can be conve- 
niently performed in hyperspherical coordinates. For this 
purpose we introduce the function £l( k \p 2 ) |j, 



sm K * ipi 



difii sin k 1 1 ip. 



/ d(p k -i exp ~y^Jt 
Jo V i=i 



where 



and 



ti = p cos ipi sin ipj if i < k, 

3=1 



fc-1 

tk = p n sin ^' 

3=1 



(A8) 



(A9) 



(A10) 



with p — p k . The exponent can be expressed in hyper- 
spherical coordinates as 



1=1 



+ p 2 sin 2 ^i(l + sin 2 ^2(l + - ■• (l+sin 2 v?fe-i)))- (All) 

Inserting this expression into Eq. ( |A8| ) and expanding 
the exponential function into a Taylor series, we arrive at 
Eq. @ together with Eq. @. A recurrence formula for 

the coefficients A n in Eq. ( |l0| ) can be readily obtained: 
AM=B 2n+k ^( n \A?-'\ fc>3 (A12) 



where 



r ( — ) 

B, = I dip sin 3 ip = Vtt 2 2 

- 1 - V 2 / 



(j > 0). (A13) 



Starting from AX> = 2B 2n , the formulas ( |A12| ) and (g) 
allow for fast and accurate numerical determination of 
the functions QS k '{p 2 ) even for high fc. 



APPENDIX B: ASYMPTOTICS OF fi (K) (p 2 ) AND 
DIVERGENCE OF KERNELS K K (R) 

In order to analyze the divergence of the kernels K k (r), 
we must first know the asymptotic behavior of the func- 
tions £l( fc )(p 2 ) for large p. We start from the in tegral 
representation (A8) and write the exponent (All) as 



H k = P 2 + P 2 sin 2 ipi$(ip 2 , • ■ • , ifk-i), (Bl) 



i=i 



with 



$(<p 2 ,...,<Pk-l) = 1 + sin 2 V3 2 (l + sin 2 ip 3 (l + ...)). 



(B2) 



Note that $ > 1. We insert Eq. (§1) into Eq. (|AJ) and 
integrate over cpi. Assuming k > 2, the relevant integral 



is 



I = 



i • k — 2 

d(pi sin y?i e 



-p <I> sin </?i 



2 / dt/^i sin (^i e 
'o 



— p $ sin t/?i 



(B3) 



Assuming p 2 3> 1, we may write sin ~ 9^1, because the 
integrand is essentially nonzero only for (pi <C 1. From 
the same argument, we can extend the integration from 
(0,7r/2) to (0,oo), 



/«2 dipx ipi ~e 
Jo 



(/$r 



(B4) 



To finish the calculation of O, one has to integrate 
<j,(i-fc)/2 over .j-ke remaining angles ip 2 , ■ ■ ■ , <pk-i (with 
appropriate measure). We eventually find the asymp- 
totic behavior 



p 



c kP 



i-fc. 



(B5) 



The factor exp(— p 2 ) comes from the first p 2 in Eq. flBl| ). 
Taking into account that fiWfp 2 ) = 2exp(— p 2 ), we see 
that the asymptotic behavior (B5) holds for all fc > 1. 

To investigate the divergence of K k (r) at r — ► 0, we 
make use of the integral representation (^), which is 
rewritten here as 







K k {r) 



T k/2 



x e {k+1)p exp 



-(e» -1) 



(B6) 



For small r, the dominant contribution to this integral 
comes fr om l arge p. We can replace f2 by the asymptotic 
formula (B5) and absorb the prefactors into C k , 



K k {r)^C k r k / dpe kp exp 
Jo 

Change of the variable according to 
t = (e p2 - l)r 2 

yields 



-{e p -l)r 2 



(B7) 



(B8) 



K k (r) 



1 

2^k 



dt In 



t + r z 



-1/2 



(t + r 2 f 



(B9) 



The integration region of flB9|) can be divided in two 



parts, t < R 2 and t> i? , with R i 



For small r, the 



dominant contribution stems from the latter part where 
the approximation 



In- 



f- 



ill- 



can be used, and we find that 



1 



(—9 In r^ 1 / 2 r°° 

K k (r) « / dt(t + r^e-K (Bll) 



(B12) 



2r k 

The integral is finite for r— > 0, and thus 

Kk{r) ~ F^iTF 1 ' 

The logarithm singularity in the denominator is very 
weak in comparison to the polynomial one in the nu- 
merator. Only the polynomial divergence is relevant in 
Eq. ([!]), because it determines whether the integral is 
convergent or not. Thus we ne ed not consider the loga- 
rithmic part, so that Eq. (B12) simplifies to 



K k (r) 
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